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Abstract. The space of phase-parameters (sky-position, frequency, spindowns) 
of a coherent matched-filtering search for continuous gravitational waves from 
isolated neutron stars shows strong global correlations ("circles in the sky"). 
In the local limit this can be analysed in terms of a parameter-space metric, 
but the global properties are less well studied. In this work we report on 
our recent progress in understanding these global correlations analytically for 
short to intermediate (less than a month, say) observation times and neglecting 
spindowns. The location of these correlation-circles in parameter-space is found 
to be determined mostly by the orbital velocity of the earth, while the spin- 
motion of the detector and the antenna-patterns only contribute significantly to 
the amplitude of the detection statistic along these circles. 
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1. Introduction 

Continuous gravitational waves from, for example, rotating neutron stars with non- 
axisymmctric deformations such as mountains or oscillation-modes, are one of the 
primary targets of current-generation detectors of both interfcrometric- (e.g. see [2, 3]) 
as well as bar-detector design [9], 

In this paper we restrict ourselves to gravitational waves emitted from isolated 
neutron stars with negligible proper motion, for which the signal can be assumed to 
be a nearly monochromatic sinusoidal (with slowly decreasing frequency) in the solar- 
system barycenter frame (SSB). The corresponding signal received at the detector will 
be Doppler-modulated by the spin of the earth and its orbital motion around the sun. 
The phase of the received signal therefore depends not only on its intrinsic frequency 
/ and spindowns = d k f/dt k , but also on the sky-position (denoted by a and 5 
for right ascension and declination) of the source and on the detector location. In 
addition there is a time-dependent amplitude modulation of the signal depending on 
its polarisation angle ip and amplitudes h + and h x . However, as shown in [5], one can 
eliminate these "amplitude parameters" together with the initial phase of the signal 
$o by analytically maximising the detection statistic over these parameters. The 
resulting reduced parameter-space therefore only consists of the "phase parameters" 
V = {a,5, f, f( k '}, and the corresponding partially maximised detection statistic is 
usually referred to as the "^"-statistic" . 
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Given a stretch of data x(t) from the detector, an untargeted search (as opposed 
to a search for known pulsars) consists of calculating the detection statistic T over 
the parameter space (by sampling individual points) and determining the location 
of "candidates" that cross a predetermined threshold. For a given signal V s = 
{a s , S s , f s , / s }, the detection statistic T in a search-point V generally depends on 
the parameter-mismatch AV = V s — V, i.e., T = !F(V S , AV) (strictly speaking, T 
also depends on the amplitude parameters of the signal). Neglecting noise, one can 
expect T to be maximal for the perfectly matched search-point AV = 0. In the local 
neighbourhood of this point, the loss AT = T(V s ,0)—T(V s , AV) in detection statistic 
can be quantified in terms of a parameter-space metric as discussed in [8, 4]. Namely, 
neglecting effects of the amplitude-parameters of the signal, this metric can be defined 
as AT(V S , AV)/T(V S , 0) = £^ gij (V s ) AV 1 AV j + 0(AV 3 ). 

In the case of isolated neutron stars, we find that this metric is highly anisotropic 
in terms of the above parameter-space variables (in particular with respect to sky- 
position), as will also be seen in section 6. Furthermore, the global behaviour of 
T differs dramatically from the local quadratic decrease as a function of parameter- 
mismatch AV. This can be seen, for example, in figure 2 for an all-sky search at 
fixed frequency of an injected signal. For each target search-frequency / within the 
Doppler-window of the signal frequency / s , one finds a very thin circular band on the 
sky along which T is of comparable magnitude to its maximum, while it drops sharply 
to zero in the directions orthogonal to this "circle" . In the local limit this feature has 
already been pointed out in [7] . 

In this work we present an approximate analytic description of these "circles in 
the sky" , tracing their physical origin to the Dopplcr-shift due to the orbital motion 
around the sun. The Doppler-shift induced by the spin-motion of the earth as well 
as the amplitude modulations due to the rotating antenna-pattern are found to give 
only very small corrections to the structure of these patterns in parameter-space (i.e., 
the location of the circles) , while they do contribute significantly to the amplitude of 
the ^-"-statistic along these circles. 

2. Matched filtering of continuous signals 

The gravitational- wave strain h(t) at the detector can be written in the form 

h(t) = F+(t) h+(t) + F x (t) h x {t) , (1) 

where F +tX are the antenna-pattern functions, and h + x are the two polarisation 
components of the gravitational wave. For a continuous pulsar-signal we can write 

h+(t) =A + sin($(t)+$ ) , h x {t) = A x cos ($(<) + * ) • (2) 
The phase $(£) of the signal is assumed to be of the form 

S»(t; / <*> , n) = 2n £ y^yy r fe+1 (t, ft) , (3) 

where s is the number of spindown parameters, and r is the arrival time 
in the solar-system barycenter (SSB) of a wave-front from direction n = 
(cos 5 cos a, cos S sin a, sin 5) arriving at the detector at time t, i.e., 

/ f(t) ■ ft 
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with r(t) denoting the vector from the SSB to the detector. As shown in [5], the full 
phase-model (3) can be simplified by an approximated model without significant loss 
in detection statistic, namely 



$(t) « 2tt 



y / (fc) 

^ (k + 1) 



fe+i 



.k=o 



k=0 



(5) 



which will be dominated by the first term, i.e., $(t) s» 2tv ft. In the case of stationary 
Gaussian noise, the likelihood ratio A (which is the optimal detection statistic in the 
sense of Neyman-Pearson) is given by 



T 

lnA = — 

Oh 



{x\\h)-\{h\\h) 



(6) 



where Sh is the power spectral density of the noise, and the scalar product (x\\y) is 
defined as 

(x\\y) = ~ f T x(t)y(t)dt. (7) 



3. Simplified matched-filtering statistic: neglecting amplitude 
modulations 



In the following we make the simplifying assumption that we can neglect the 
amplitude-modulation caused by the antenna-pattern, so we assume F + ^ x « const., 
and using (1) and (2) this results in the simplified strain model 

h(t) =Ai cos$(t)+A 2 sin$(i). (8) 

Maximising the log- likelihood function (6) over the two unknown amplitudes Ai^, we 
obtain 

^ A Mi = ^-\(x\\e- i *)\ 2 = ^-\X\ 2 , (9) 

defining the matched-filtering amplitude X. In the following we assume the data 
contains only a unit-amplitude pulsar-signal of the form 



x(t) = n[s(t)], with s(t)=cxp i^^/fUs) 



characterised by the signal phase-parameters fs k ^ and n B . 
amplitude X, defined in (9), can therefore be expressed as 



(10) 

The matched-filtering 



X = die 



(11) 



where A$ is the phase-difference A$(i) = <1> S (t) — <&(t), and where we have used the 
fact that J exp[i27r(/ + f s )t]dt w 0. We see from (11) that \X\ will have a global 
maximum of |A| = 1 if the phase-parameters of the template are perfectly matched 
to the signal, i.e., A$ = 0, while \X\ will decrease very rapidly for increasing phase- 
mismatches. 

The central aim of the following investigation is to determine the regions in 
parameter-space for which the detection statistic \X\ is of order unity for a given 
signal {n s , /J fc) }. 
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4. Analytic approximation for \X\ 2 

In order to simplify the following discussion, we restrict ourselves in this study to 
pulsar-signals with negligible spindown over the observation period T, so we assume 
/( fe ) k for k > 0. Using (5), we can therefore write the phase- mismatch A$(i) as 



c 



(12) 



A$(t;ra a ,/ s ;n,/) = 27r 
where wc defined 

A/ = /„-/, and A = f s H s -fH. (13) 

We can decompose the motion of the detector r(i) into its orbital- and spin-component, 
namely, f(t) = r or t,(t) + F sp i n (t), where r or b is the vector from the SSB to the centre 
of the earth, and F sp ; n is the vector from there to the detector. Correspondingly we 
have 

A$(t) = 2nAft + A$ OTh (t) + A$ spin (i) , (14) 

where 

27T 27T 

A$ orb = — r orb (t) ■ A , and A$ spin = — f spin (t) -A. (15) 
4-1- Spin motion of the earth 

Except for very short observation times T<Cl day, which are not relevant for searches 
for continuous waves, the phase modulation $ sp in due to the spin-motion of the earth 
is oscillatory, and can therefore not be treated using a Taylor-expansion. However, 
following an approach first used in [6, 10], it is more advantageous to use the Jacobi- 
Anger expansion, which allows one to expand an oscillatory exponent in terms of the 
Besscl- functions J n (z), namely 



e iz cos C _ >^ 



i n J n (z)e"<. (16) 



n— — oo 

In order to be able to use this, wc write the detector motion r sp i n (i) as 

r sp in(t) = sin A R® z + cos A i?® m(t) , (17) 

where A is the latitude of the detector, is the radius of the earth, z is the unit- 
vector along the rotation axis and w{t) is the time-dependent unit-vector orthogonal 
to z, pointing from the rotation axis to the detector. Using this decomposition, we 
can write the phase-mismatch due to the spin-motion of the earth as 

A$ spin (£) = A$, + K A cos <p(t) , (18) 

where wc defined 

A$ z = 2tt— sinAA z , K A = 2n— cos A A x , (19) 
c c 

and 

tp(t) = ip A + fljpin t , with <p A = z(A±,m(0)\ , (20) 
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and where A z and A± are the parallel and orthogonal projections of A with respect 
to the rotation axis z. Using this decomposition, we can apply the Jacobi-Angcr 
expansion (16) as follows: 



iA* Hpin (t) 



i A<I>~ 



J2 i n e mVA J n (K A )e tnn "^ t . (21) 

n— — oo 

Substituting this together with (14) into the matched-filtering amplitude (11), we 
obtain 



X 



e 



J2 i n e in ^J n (K A )Y n , 



n— — oo 



where 



Y n = - J exp [i (2ttA/ t + A$ orb (i) + n fi spin t)] dt . 



(22) 
(23) 



4-2. Orbital motion around the sun 

Contrary to the spin-motion, the orbital motion induces a slowly varying and non- 
oscillatory phase-correction for observation-times T < 1 year, so we can expand the 
orbital phase- mismatch A<J> or b(i) in terms of the small quantity e = Q ibT <C 1, which 
leads to 



with the coefficients 
ai 



D iA* orb (0) 



1 



exp 



i ai 



t + a kt k I 

k=2 ) 



dt . 



— = Af + /3 -A + nf spin , 

Z7T 



(fc>l): 



(24) 



(25) 



where we defined / sp in = ^spin/27r, and where /?o = Vq/c = r O rb(0)/c denotes the 
orbital velocity at the start of the observation t = 0. 

Note that the integrals Y n will be strongly peaked as a function of the index n, 
namely by considering 



dt 



T 



(26) 



we see that Y n will have a maximum for the index n* closest to the point a\ = 0, i.e., 

~Af + o -A~ 



= —round 



/spi] 



(27) 



As a first approximation we therefore only keep the dominant term n* in the Jacobi- 
Angcr expansion, i.e., we set Y n ss Y m and so we obtain for the matched-filtering 
amplitude 



X 



2 Jn* (-^Ca) ^n* 



(28) 



We are interested in the regions of parameter-space where \X\ ~ 0(1). In these regions 
we can assume the exponent in (24) to be small, and therefore, after Taylor-expanding 
and integrating (24), we can obtain the approximate expression 



\X\ 



where d\ = ai(n*). 



1 _ £l t 2 - — 2 rjl3 1 
12 6 



■T 3 + 0(T 4 ) 



= \Xi\ 



(29) 
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5. The 'maximum structure' of \X\ 

On can sec from its definition (11) that \X\ has a global maximum of \X\ = 1 in 
the perfectly matched case (/ = / s and ft = n s ), while it will be nearly zero almost 
everywhere else, because generally A$ will not be small. Guided by the observed 
structure of the detection statistic (see section 6, for example), we focus on the points 
in parameter space where \X\ (or T) has a local maximum in target-frequency / for 
a fixed sky-position ft, i.e., 

In order to understand the qualitative features of these local maxima, we now also 
neglect the spin motion of the earth (setting n* = and = 0), and using (29) and 
(25) we obtain the following condition for local maxima (in the sense of (30)): 

a\a\ + a-i Ta[ + ai Ta' 2 = , (31) 

where we write a' k = dfak(f,ft). We can solve this order by order in T (strictly 
speaking, in the small parameter e = f2 rt>T), and to zcroth order we find 

O{T ) : a x a[ = => A/ + (3 ■ A = , (32) 

as a[ = 0(1). For the next order we therefore solve (31) with the ansatz a\ = XT 
(with A = 0(1)), so we find 

0{T) : ai +a 2 T = =^ Af + 0i-A = O, (33) 

where /3i = (3q + (3 -j is the first order approximation of the orbital velocity at the 
midpoint of the observation time. This scheme can be extended to arbitrary order 
in T, but for the qualitative discussion here we restrict ourselves to first order. The 
condition (33) can be expressed more explicitly as 

/ (l + fk ■ n) = f s (l + fk ■ rt s ) . (34) 

For a given signal {n s , / s }, this equation describes the (approximate) location of the 
local maxima of \X\ as a function of / for fixed ft. Conversely this can be regarded 
at fixed frequency /, describing the sky- locations ft in which the maximum occurs at 
this particular frequency /. From the above equation one can easily see that these 
points ft describe a circle in the sky, which has its centre in j3\, and an opening angle 
determined by the signal- location and the target frequency. In the case / = / s , we 
see that ft = ft s is one of the possible solutions of (34), and so the circle will pass 
through the signal- location in the sky. If / < / s , then we must have (3\ ■ ft > f3\ • H s , 
i.e., the opening angle of the circle has to be larger than in the case of a perfectly 
matched target-frequency, and vice- versa for the case of / > / s . For a 10- hour 
integration starting at GPS-time £o = 732490940 s, the orbital velocity at the midpoint 
is (3i = [0.36, —9.13, —3.96] x 10 -5 , and the corresponding family of circles is shown 
in figure 1 . This equation also determines the upper limit on the frequency-mismatch 
A/ for which such a local maximum is possible, namely, by rewriting (34) as 



A/ 



ft • An 



1 + A • n 



<-^- = 2f3 1+ 0({3t), (35) 
1 - Pi 
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-1.5 -_ 



1 2 3 4 5 6 

Right Ascension [rad] 

Figure 1. Family of circles described by (34) for a T = 10 hours observation 
starting at GPS-time to = 732490940 s. The direction of 0\ is marked by 'x', 
while '+' indicates the direction —f)\. 

we see that the maximal relative frequency-mismatch is bounded by twice the Dopplcr- 
shift due to the orbital velocity. Using the value V OI b w 3 x 10 4 m s _1 , this corresponds 
to a "Doppler-window" of Af/f ~ ±2 x 1CT 4 . 

We note another interesting property of the approximate expression (29) for the 
matched-filtering amplitude \X\: namely, from the above discussion we know that 
the local maxima of \X\ occur on or close to the circles described by (34), which are 
characterised by a\ k. 0. Therefore we see from (29) that the amplitude of |Xi| along 
these circles will be dominated by the factor J„*(i^^), where Ka is proportional to 
A±, the projection of A into the equatorial plane. This factor is therefore symmetrical 
in the target sky-position ft with respect to the equatorial plane, a feature which can 
indeed be seen in the plot of (29) shown in figure 4(b). Interestingly, this symmetry can 
still be present in the fully numerical .7-"- statistic, as seen, for example, in figure 3(b), 
while the effects of spin-motion and antenna-patterns can also mask this feature 
depending on the parameters, as seen in figure 3(a). 

6. Comparison to fully numerical results for the JF-statistic 

We now turn to comparing these analytic results to the actual parameter-space 
structure of the fully numerical ^-statistic including all effects of amplitude- 
modulation and using an ephemeris-based description of the detector-motion. As 
an example, we consider a 10-hour observation time starting at GPS-time to = 
732490940 s, and we use a fake pulsar-signal at / s = 100 Hz and sky-position 
(a s , 5 S ) = (2,1), for two different choices of amplitude-parameters shown in 
table 1. In the following the detector- location is always chosen to be the 
LIGO Livingston Observatory (LLO), and we only consider pure signals without noise. 
The fake pulsar-signal is generated using makef akedata_v4, and the J-"-statistic [5] is 
calculated using ComputeFStatisticl 

X These codes are found in LALApps[l], under src/pulsar/Injections and src/pulsar/FDS_isolated. 
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Table 1. Amplitude- and phase- parameters of two fake pulsar-signals. 



Model 


/ [Hz] 


(a, S) [rad] 


A+ 


A x 


* 


<1>() 


I 


100 


(2.0, 1.0) 


1.0 


0.5 


1.0 


2.0 


11 


100 


(2.0, 1.0) 


1.0 


0.0 


0.0 


0.0 



We observe the following structure of the ^-statistic over the parameter space 
{a, 6, /}: searching over the sky at a perfectly matched target-frequency / = / s , the J-- 
statistic has a maximum at the signal sky-position (a s , 5 B ), as expected, but remains of 
the same order on a very narrow, complete circle over the sky (obviously including the 
signal-position), while it is practically zero everywhere else. For mismatched target- 
frequencies / (within the Doppler- window of / s ), there is a different but concentric 
circle of J-- values comparable to the maximum, with a larger radius for f < f s , and 
a smaller radius for / > / s , as shown in figure 2. This agrees qualitatively and 
quantitatively surprisingly well with the circle-structure predicted solely based on the 
effect of the orbital Dopplcr-shift, i.e., (34) and figure 1. 




-1.5 b | I I I L_ 

1 2 3 4 5 6 



Right ascension [rad] 

Figure 2. Plot of the jF-statistic over the sky superposed for different target 
frequencies /: (1) / = 100.002 Hz, (2) / = 100.001 Hz, (3) f = f B = 100 Hz, 
(4) / = 99.997 Hz and (5) / = 99.99 Hz. The parameters of the injected signal 
correspond to model I in table 1, and the small square indicates the sky-position 
of the signal. 

In the next step we run an all-sky search over the whole Doppler-window 
/ € (1 ± 2 x 10~ 4 )/ s , corresponding to a dense superposition of these circles for 
the different target-frequencies. The resulting ^"-statistic projected on the sky is 
shown in figure 3(a) and (b) for pulsar-signals I and II, respectively. From these 
projected plots we see more clearly that the variation of the amplitude of the T- 
statistic maxima is more sensitive to the polarisation-parameters of the signal, due to 
effects of the antenna-pattern. For example, figure 3(b) shows some of the equatorial 
symmetry characteristic for the approximate expression (29), as discussed in section 5, 
in particular we see a mirror-image in (2, —1) of the absolute maximum at the signal- 
position (2, 1). For signal I, on the other hand, the projected T is more asymmetrical 




Right Ascension [rati| R|ght Ascension [rad] 



Figure 3. Projected jF-statistic over the sky in the frequency range / S 
(1 ± 2 X 10 — 4 ) / a for (a) pulsar-signal I and (b) signal II as denned in table 1. 

and there is no mirror- image of the maximum in the southern hemisphere. 

Summarising these observations we can state that the locations of the local 
maxima of T are dominantly determined by the orbital velocity of the earth, and 
are described very well§ by the circles (34), independently of detector- location and 
amplitude-parameters of the signal. The spin-motion of the detector and amplitude- 
modulation due to the rotating antenna-pattern, on the other hand, are responsible 
for the variation of the amplitude of T along these circles, as seen in figure 3. 

Finally, we plot the approximate detection amplitude |^fi| 2 given by (29) over the 
sky for different target frequencies, as shown in figure 4(a), which shows some striking 
similarities to the 'exact' result in figure 2. In particular, it exhibits the feature 
of vanishing very rapidly outside the narrow "circles" described approximately by 
(34). The amplitude of |^i[ 2 on the circles, however, is not a good description of the 
corresponding ^"-statistic amplitude, as |^i| 2 is seen to be much more strongly peaked 
around the signal- location and decreases rather rapidly to about 0.1 on most of the 
circles (see figure 4(b)), while the ^-"-statistic has a similar amplitude to the maximum 
over most of the circles. Note also that figure 4(b) illustrates the equatorial symmetry 
of \Xx | 2 as discussed at the end of section 5. 

7. Discussion 

In a practical search for continuous gravitational waves we will usually be interested in 
the highest values ("candidates") of T over a chosen volume of the parameter-space. 
From the results presented here we expect that a real signal would generate candidates 
nearly everywhere on the sky within the Doppler-window ±2 x 10~ 4 of the signal- 
frequency / s . All these candidates, however, will satisfy (34) with an (approximately) 
identical value on the right-hand side, and so the quantity n = /(l + j3i • n) is an 
(approximate) invariant of candidates caused by the same signal. This could be a 
very useful criterion for classifying equivalent candidates, and also for determining 

§ A preliminary numerical comparison of the frequency fjr of the maximum of T for a fixed sky- 
position ft and the theoretical prediction f th from (34) yields the bound \fjr — /thl//^ < 10~ 6 over 
the whole sky. 
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Right ascension [rad] 



Figure 4. Plot of the approximate matched-filtering amplitude |Xi| 2 over the 
sky for (a) different target frequencies (1) / = 100.002 Hz, (2) / = 100.001 Hz, 
(3) / = /«. = 100 Hz, (4) / = 99.997 Hz and (5) / = 99.99 Hz, and (b) over the 
whole Doppler-window / £ (1 ± 2 X 10 _4 )/ s projected on the sky. The small 
square indicates the sky-position of the signal (model I in table 1). 



coincident candidates from different detectors. One could further use this to construct 
coincidence-tests of candidates from different observation times, as (3\{t) is a known 
function of time. In order for these methods to be practically useful, however, the 
errors of the different approximations need to be quantified, and we need to test the 
robustness of the suggested coincidence criteria in the presence of noise. 
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